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In the present work, we revisit the highly active research area of inhomogeneously nonlinear 
defocusing media and consider the existence, spectral stability and nonlinear dynamics of bright 
solitary waves in them. We use the anti-continuum limit of vanishing coupling as the starting point 
of our analysis, enabling in this way a systematic characterization of the branches of solutions. Our 
stability findings and bifurcation characteristics reveal the enhanced robustness and wider existence 
intervals of solutions with a broader support, culminating in the “extended” solution in which all 
sites are excited. Our eigenvalue predictions are corroborated by numerical linear stability analysis. 
Finally, the dynamics also reveal a tendency of the solution profiles to broaden, in line with the 
above findings. These results pave the way for further explorations of such states in discrete systems, 
including in higher dimensional settings. 


I. INTRODUCTION & BACKGROUND 

Intrinsic localized modes have received considerable attention during the past two and a half decades since their 
theoretical inception; see e.g. During this time frame, a wide range of experimental contexts has been shown 

to support their existence and robust role in the systems’ dynamics. These range from arrays of nonlinear-optical 
waveguides Q to Bose-Einstein condensates (BECs) in periodic potentials Q, and from micromechanical cantilever 
arrays @ to Josephson-junction ladders Additional areas of application include (but are not limited to) granular 
crystals of beads interacting through Hertzian contacts Q, layered antiferromagnetic crystals Q, as well as halide- 
bridged transition metal complexes 0, and dynamical models of the DNA double strand [l^. This broad span of 
activities has been summarized in [^, although admittedly applications keep being added every year; see e.g. as 
recent examples BUI- 

A model that has contributed considerably to advancing our understanding of such lattice nonlinear dynamical 
systems and their coherent structures is the so-called discrete nonlinear Schrodinger (DNLS) equation [l^ . Its appar¬ 
ent simplicity as a prototypical system incorporating the interplay of nonlineariW and a discrete form of dispersion, 
but also its relevance as a suitable approximation of optical waveguide systems [3, ll^ s-nd atomic systems in optical 
lattices 0 have, undoubtedly, contributed to the popularity of the model. Another key feature is its ability to capture 
numerous linear and nonlinear features observed in experiments such as discrete diffraction [l^ and its manage¬ 
ment , discrete solitons [11,111 and vortices B, Bl , Talbot revivals [2l| , and VT -symmetry and its breaking Bl , 
among many others. 

On the other hand, a theme that has met with growing interest in nonlinear Schrodinger settings recently is that 
of spatially modulated nonlinearities; see e.g. for a review geared specifically towards periodic modulations the work 
of [^. Within that theme, a focal point has been the study of the potential bright solitons in the context of defocusing 
nonlinearity, as supported by spatial modulations. This possibility introduced for solitons and vortices in I24l - l26l| was 
also extended to the setting of spatially inhomogeneous nonlinear losses in , quintic nonlinearities in l28ll . domain 
walls in (^|, F ermi and Bose gases in 1^. dipolar Bose-Einstein condensates in [3l|, nonlocal media [32j| . discrete 
systems in 33| and even in 3d media in [34l| . 

In the present work, we combine the two above themes. In particular, we explore the defocusing DNLS equation 
in the presence of a growing nonlinearity. Contrary to the case of [3^ . here the nonlinearity does not present an 
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exponential growth, but rather a power law in our concrete implementation. Nevertheless, this is not the focus of our 
contribution. Instead, we keep the analysis as general as possible, considering an arbitrary profile in the nonlinearity 
(given by g{n)) whenever possible. Our aim is to start from the well-established anti-continuum limit of and 
following the existence and stability considerations of [s^ to provide a systematic view of the possible excitations 
in the discrete system in the form of bright solitary waves. We examine states with 1, 2 and 3 sites, as well as the 
“extended” state where all the sites are excited. We reveal the stability of the different states, and also explore how 
approaching the continuum limit, more and more extended states are favored (while “narrower” states disappear in 
suitable bifurcations), leaving as most suitable state in the limit the extended state mentioned above. We also provide 
a comparison of the stability properties with the corresponding homogeneous limit (where the nonlinearity strength 
is equal for all sites), illustrating that the inhomogeneous nonlinear profile effectively promotes the instability of the 
few-site localized states considered. 

Our presentation will be structured as follows. In section 2, we provide the theoretical background, analyzing 
the existence and stability properties of different states. In section 3, we provide numerical existence/linear stability 
results that corroborate the analysis, as well as direct numerical simulations following the dynamics of unstable states. 
Finally, in section 4, we summarize our findings and present a number of directions for future consideration. 


II. THEORETICAL ANALYSIS 


The model that we will consider will be of the form: 


iiln = -e {Un-l + Un+i - 2u„) + g{n)\unfun. (1) 

Here, we have in mind a defocusing nonlinearity, as is e.g. the case in LiNbOs [svl - fs^ . used previously to demonstrate 
a number of features including dark and dark-bright solitary waves. However, the waveguides in the present setting 
are effectively “tailored” to have distinct Kerr response, forming the profile associated with g{n) > 0. 

We will seek standing waves in the customary form Un = {t here plays the role of the -spatial- evolution 

variable and A is assumed to be positive), hence will satisfy: 

AVn = -eA2Vn+gin)\Vn\'^Vn. ( 2 ) 

Here A 2 plays the role of the discrete Laplacian. At the anti-continuum (AC) limit of e —>■ 0, the only solutions 
are = 0 and = \JAjg{n)e''^'^. Enforcing the latter for every site and provided that g(ji) grows indefinitely leads 
to a decaying pulse (the extended solution considered herein). While we will briefly touch upon this waveform, our 
emphasis will be (similarly to [^ 1 on few-site excitations. 

Using Eq. ([2|) multiplied by the conjugate u* and subtracting from it the conjugate of Eq. ([2]) multiplied by we 
obtain a solvability condition which is the same as in the standard DNLS case, namely: 


0 = • • • = - u>„+i = • • • = 0, 


( 3 ) 


due to our implicit assumption of |u„| —>■ 0, as n —>■ 00. Using the AC limit solution of = -\/A/g(n)e®®”, this yields 
that for adjacent excited sites the condition sin(0„+i — On) = 0 should hold allowing relative phases only of 0 or tt for 
such sites. 

We now explore the corresponding linearization problem, using the ansatz 


Un = e 


—iAt 


[pne^* +qn 


,A t ' 


( 4 ) 


(where * denotes complex conjugate) and deriving the equations to 0{6) for (a„, bn) such that Pn = Un + ibn and 
Qn = Un — ibn, given the complex nature of the perturbations to the solution Vn (36|. Notice that hereafter, we will 
restrict ourselves (without loss of generality for the one-dimensional setting) to real solutions, assuming = 0 or tt. 
Here we are effectively using the gauge invariance of the DNLS to fix one of the excited sites’ phase to 0 (or tt) and 
the solvability condition above to obtain that all other excited site phases should then also be 0 or tt [T^. Then, the 
resulting eigenvalue problem reads: 


A 


(In 

br) 




( 5 ) 


Here C-bn = —eA 2 bn — Abn + g{n)Vnbn, while £+a„ = —eA 2 a„ — Aa„ + 3(7(n)u^a„. Combining the two linearization 
equations, we obtain 


X^bn = -C+C-bn ^ X^C^bn = -C-bn- 


( 6 ) 
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Notice that in the vicinity of e —^ 0, £+ becomes a multiplicative operator with positive entries, hence is invertible. 
Now, forming the inner product with 6 „, we obtain 

^2 _ {bn, ^-bn) 

But again, near the AC limit —>■ A/g{n), leading to £+ —>■ 2A and hence —>■ (2A)“^. Therefore, the stability 

will critically hinge on the eigenvalues of 

£- can be directly seen by the considerations above to vanish at the AC limit for the excited sites. For the non- 
excited sites., = C- = —A, yielding A = ±Ai. Hence, at the AC limit, except for the N excited sites corresponding 
to 0 eigenvalues, all other eigenvalues will be degenerate at ±Ai and as e becomes nonzero will form the continuous 
spectrum [A — 4e,A]; hereafter and without loss of generality we will set A = 1. However, the key for stability 
considerations regards the TV — 1 eigenvalue pairs bifurcating from the origin in the case of N excited sites (one pair 
will stay at A = 0 due to the gauge invariance of the model). To determine these eigenvalues, it is critical to evaluate 
the N X N reduction of the operator £_ so as to obtain the eigenvalues from Eq. 0 . To do so, we follow a similar 
approach as in [s^ expanding Vn = Vn'^ + evn'^ + ... and computing the leading order correction as: 

_ 1 f cos( 0 n+i - 0 „) , cos( 0 n-l - 0 n)] i 0 ^ 

" “n 

when for the n-th site both of its neighbors are excited; when only one neighbor is excited, then only the corresponding 
term is present in Eq. Notice that in this expression and hereafter for simplicity (and without loss of generality), 
we will set A = 1. Using this expression the diagonal elements of the N x N matrix arising in the numerator of Eq. ([7]) 
M = {b,C-b) are found to be: 


Mn,n = \/ gin) 


C0s(6ln+1 - On) 

\/9{n+l) 


C0S(6>„_1 - 0n) ^ 

Voin - 1 ) / 


(9) 


(again, if both neighbors are excited). On the other hand, the off diagonal contributions remain the same as in [s^ . 
namely 




( 10 ) 


Once the eigenvalues 7 of A4 are calculated, then the eigenvalues of the full problem bifurcating from 0 are given 
as A = ±y/—2e’~f. Let us give some explicit examples. In the case of iV = 2 excited sites, the relevant matrix 


M = 



C0S(6'„+1 - 0n). 


( 11 ) 


This leads to 7 = 0 and 7 = ) cos(0„+i — 0„). It is particularly interesting to note that in this 

setting (+ \j with the latter value being the homogeneous limit case of g{n) = 1 (i.e., of all sites 

bearing an equal nonlinearity prefactor). This effectively implies that the inhomogeneous solution will always be more 
prone to instability. In the out-of-phase case of cos(0„+i — 0n) = —I, this will be because of a real eigenvalue pair 
which is larger in magnitude in the inhomogeneous case. On the other hand, in the in-phase case of cos(0„+i —0n) = Ij 
the eigenvalue pair will be imaginary (again larger in magnitude for the inhomogeneous case) and will start growing 
along the imaginary axis as e is increased. This, in turn, given (as in the case of [ 3 ^ : see the relevant discussion 
therein) the negative signature of the relevant eigenvalue, will eventually lead to an instability as e increases, upon 
the collision of this eigenvalue with the continuous spectrum; see also below the detailed discussion associated with 
Fig. m Based on the above discussion, we expect the inhomogeneous case to be more prone to instability than its 
homogeneous counterpart. 

In the case of TV = 3 excited sites, the resulting 3x3 reduced matrix is of the form: 


M = 


^ COS(0„_1 0n) 

- cos( 0 „_i - On) 

0 


0 ^ 

- cos( 6 '„+i - On) (12) 
^B^^COs{0n+l - On) J 


V 
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In this case too, one can find explicitly the eigenvalues of the matrix, although the resulting expression is 
far more cumbersome. More specifically, setting a = b = iii addition to 0, the 

other two resulting eigenvalues 7 in this case are: 7 = ( 2 a 6 )*^“^)((l + a^)cos( 0 „_i — 9n) + acos( 0 „+i — 0„) ± 
A/-4a6(&2 + a?{l + b"^)) cos( 6 >„_i - 6 <„) cos( 6 >„+i - 6 »„) + ((1 + a^)bcos{9n-i - 6n) + a(l + b'^) cos( 6 »„+i - 6»„))2). In 
this case, it is less straightforward to provide a general statement about the comparison to the homogeneous state of 
g{n) = 1. Nevertheless, we would like to note that in the case considered here also, the magnitudes of the eigenvalues 
are larger in comparison with the homogeneous case of g{n) = 1 Vn (for the corresponding setting of three-site- 
excitations within the latter). This, in turn, provides a stronger (larger growth rate) instability -again in comparison 
with the homogeneous g(n) = I Vn case- for configurations with at least an out-of-phase pair of adjacent sites and an 
instability arising for smaller values of e in the case of (all) excited sites bearing an in-phase structure with respect 
to their neighbors. 

Lastly, although we give no quantitative information about that case, it is relevant to add a brief note regarding 
the extended excitation. In the latter case, it is important to point out that all eigenvalues are at 0 in the AC limit. 
Hence, the size of the matrix A4 in this case is comparable to the domain size and hence it is less straightforward 
to characterize the relevant eigenvalues. On the other hand, it is especially relevant to report that since the pulse 
is decaying in the case of a potential growing at infinity, the corresponding effective potential will be unbounded as 
n —>■ 00 in that case. This, in turn, leads the eigenvalues to bifurcate from the origin of the spectral plane with A = 0 
giving rise to a point (rather than continuous) spectrum. We now turn to the numerical examination of the relevant 
findings. 


TABLE I: Existence intervals for the solutions considered herein. The first column labels the branches, while the second 
provides their profile form near the Anti-Continuum limit. The third column provides the end point of their termination (for 
branch G for the coupling values considered herein, no such end point was identified, hence the N/A symbolism). Finally, the 
fourth column illustrates the fate of the branches i.e., the nature of the bifurcation and with which branch they collide. It 
should be highlighted that as the bifurcation is approached, the shapes of the two (or more) colliding branches become fairly 
similar i.e., the deformation of branch B bears resemblance to branch E, and so on. In the case of branch D the collision occurs 
with the more extended branch (0,.. ., 0, ^~\f^' 0,... 0). 


Label 


Structure 


Terminal Point 


A 

B 

C 

D 

E 

F 

G 


(0, . . . , 0, 7^, 0, . . . 0) 

(0,..., 0, y^7(0), y^ 7 {T) > 0) • • • 0) 

(0,..., 0, y^ 7 ( 0 ), —\ J 7 P 7 ,0,... 0) 

(0, ...,0, ^ y^, y^^,0, ...0) 

(0,... ,0,.. .0) 
(0, ...,o,y^.-y^^^o, ...0) 

_Un = \/l/ff(w)_ 


e = 0.095 
e = 0.091 
e = 0.095 
e = 0.121 
e = 0.091 
e = 0.095 
N/A 


Endpoint Bifurcation 
Double Pitchfork with C, F 
Saddle-Genter with E 
Double Pitchfork with A, F 
Saddle-Genter 
Saddle-Center with B 
Double Pitchfork with A, G 
_ NVA _ 


III. NUMERICAL COMPUTATIONS 

For our concrete numerical example, we will use a power law growth of the nonlinear prefactor in the form g{n) = 
1 + lOn^. This choice is made purely for purposes of illustration, as the above general theory, in principle, enables the 
computation of the relevant states and their linearization eigenvalues for arbitrary forms of g(n). Arguably, the most 

fundamental branch of numerical solutions is the one with a single site excitation (0,..., 0, 0,... 0) at the AC 

limit (of e = 0). We will only touch upon this branch of solutions briefly at present and return to it, as we consider 
the bifurcations of more complex branches of solutions. In the homogeneous case where g{n) is constant, this branch 
would persist to large values of the coupling parameter forming the discrete analogue of the gap solitary wave in this 
model. Here, however, this is no longer true. Our computations show that this branch terminates around e = 0.095, 
by colliding with other solution branches as illustrated below. Interestingly, at the stability level, this branch is stable 
throughout its interval of existence. As this termination limit is approached, two eigenvalue pairs that bifurcate off 
of the continuous spectrum for e > 0.045 approach the spectral plane origin, hitting it at the critical point, a point 
indicative of the complex bifurcation scenario that will be further elaborated below. 
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We now turn to two-site solutions. We start, in particular, by considering the branch (0,..., 0, y y gW’ 0, • ■ • 0) 
that is shown in Fig. [TJ A typical solution of this sort is shown in the top right panel of the figure for e = 0.08, 
for which the corresponding example of the spectral plane (i?e(A),/m(A)) of the eigenvalues A = i?e(A) -I- ilm{\) is 
shown in the bottom right panel. The left panels in the figure illustrate the evolution of the dominant imaginary (top) 
and real (bottom) parts. We can see that as expected from the theory of the previous section and the eigenvalues 
of the matrix dm, out of the two pairs at the spectral plane origin at e = 0 (due to the two excited sites), one will 
bifurcate along the imaginary axis (blue parabolic line emanating from 0 in the top left panel) being well described 
by the theoretical prediction A ~ ±2.69y/ei (green dash-dotted line). The line of A = shows the corresponding 

eigenvalue pair for the homogeneously nonlinear case of g{n) = 1 for comparison, clearly illustrating the signihcant 
deviation of the present inhomogeneous prediction. The lower edge of the continuous spectrum is shown by the red- 
dashed line. Interestingly an eigenvalue pair bifurcating from the latter (for e > 0.04) collides with the pair stemming 
from the origin around e = 0.062 destabilizing the branch. Although this quartet briefly separates into two pairs again 
shortly thereafter, as e is (slightly) further increased the pair coming from the origin collides with the band edge of the 
continuous spectrum ensuring a quartet (oscillatory) instability for all larger values of e for which the branch exists. 

The branch appears to terminate around e = 0.091 due to its collision with (0,..., 0, — 0,... 0), 

as we will also see in what follows. 



FIG. 1: The left panels of the figure show the principal eigenvalues (top panel: imaginary part; bottom panel: real part) 
associated with the branch of solutions (0,... , 0, 0,... 0), while the right panels show a typical profile of the 

solution u„ and of its associated spectral plane (i?e(A), 7m(A)) of the eigenvalues A = Re{X) + iIm{X) for e = 0.08. In the top 
left the numerical eigenvalues stemming from the origin and from the band edge of the continuous spectrum are shown by a 
(blue) solid line. The band edge Im{X) = 1 — 4e is shown by a (red) dashed line. The theoretically predicted approximation 
of the pair bifurcating from the origin is shown by the (green) dash-dotted line, while for comparison the (lower) prediction of 
the homogeneous limit of 7m(A) = 2ye is also shown (in magenta dashed line). The bottom left illustrates that the collision 
of the pair from the origin with eigenvalues at or bifurcating from the continuous spectrum yield an instability for e > 0.062 
(see also the detailed discussion in the text). 


As the second example of a two-site excitation branch, we illustrate in Fig. [2] the out-of-phase case of 
(0,..., 0, 0,... 0). in this case the bifurcation from the origin occurs along the real (rather than 

the imaginary) axis, leading to an immediate instability of the solution. It is relevant to point out here the differences 
of this case from the corresponding focusing case; see e.g. for a relevant discussion (40j | and also for a review M- In 
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the focusing case, multi-site in-phase excitations are immediately unstable with eigenvalues bifurcating from the origin 
on the real line, while ones with all adjacent neighbors out-of-phase are linearly stable, at least close to the AC limit. 
The situation is reversed in the defocusing realm, as illustrated in this and the previous example, a pattern that will 
also be followed in the three-site excitations below. The theoretical prediction for the real pair is A ~ 2.69^, which 
we can see as a good approximation to leading order for small e, but one that progressively fails as higher orders take 
over for larger e and force the relevant pair back to the origin where it collides with an imaginary pair at e = 0.095. 
It may now be becoming clearer, both from the eigenvalue pattern and associated branch extinction collision point, 
as well as perhaps from the profile of the branch (top right of Fig. [5]) and its progressive “symmetrization” as the 
critical point is approached that this branch is involved together with the single site branch in the complex bifurcation 
further elaborated below. 



AC-limit. Here a pair bifurcating from the continuous spectrum and the real pair bifurcating from the origin (returning to it 
after an excursion along the real line) collide at the disappearance threshold of the branch e = 0.095. Again the top right and 
the bottom right panels illustrate, respectively, the (asymmetric) profile of the solution and its spectral plane for e = 0.08. 


Next we turn to three-site excitations of which (modulo permutations, similarly to [s^) we examine three examples 
in what follows. We start with the in-phase branch (0,..., 0, 0,... 0). This branch, in full 

accordance with the theory of the previous section (cf. Eq. (1121) ). has two imaginary eigenvalues bifurcating from the 
origin along the imaginary axis (while the third of the AC limit pairs of A = 0 due to the three excited sites remains 
at 0, given the phase invariance of the model). The corresponding eigenvalue pairs are theoretically predicted from 
Eq. (I12p to be A = ±2.576v^i and A = ±2.8y/ei and are found to be in very good agreement with the numerical 
findings (cf. the top left panel of Eig. [S]). Additionally, as in the two-site case, these eigenvalues are found to be 
larger than their corresponding homogeneous limit predictions, shown by the magenta dashed lines in the figure. 
This indeed implies also that the instability of the branch occurs for a smaller value of e in comparison with the 
homogeneous limit, as it arises from the collision of the two imaginary eigenvalue pairs with the continuous spectrum 
of A = 1 — 4e. In this case, given the two collisions, two oscillatory instabilities and associated quartets arise for 
e > 0.061 and 0.067, respectively. In this case, it is interesting to point out that the branch terminates for values 
larger than any of the above (as well as below) few site excited branches. More specifically, it collides with the 5-site 

branch (0 ,..., 0, — \J 0 ,... 0) at e = 0.121. This is part of a more general trend 

that we will also discuss below in the context of the extended solutions. In particular, the more extended a solution 
is, the larger the critical threshold value for its termination. 
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FIG. 3: Same as with the previous branches, but now for the solution branch with the profile 

(0,. .., 0, \/g^’ Q’ ■ ■ ■ Q) limit. The main difference here is that there are two imaginary 

eigenvalue pairs bifurcating from the origin, and two associated quartets of eigenvalues arising beyond e = 0.061 and e = 0.067, 
respectively. 


The fourth branch is the “mixed” branch with the AC limit form: (0,..., 0, 0,... 0). This 

branch has a real and an imaginary pair of eigenvalues with A = ±2.685-ye and A = ±2.685\/ez, respectively; the 
real eigenvalue pairs renders the branch generically unstable (similarly to its homogeneous counterpart) although an 
additional oscillatory instability arises for e > 0.065. We can see that the numerical imaginary eigenvalue is very 
accurately predicted by the theory. For the real one, on the other hand, we again observe the familiar feature of good 
agreement for small e, but then as e increases, higher orders take over leading to a maximal excursion along the real 
line and a return to 0 around e = 0.091 which is the termination point of the branch in a saddle-center collision with 


( 0 ,..., 


r,0,...0). 






0 ReW 


FIG. 4: Same as for the previous branches. 


but now for the “mixed” solution with profile 


(0,. .., 0, \llU)' 0,... 0) in the AC-limit. Here, one of the eigenvalue pairs moving off of the origin 

for e 0 moves along the real and one along the imaginary axis. 

Finally, from the point of view of few-site excitations, we explore the out-of-phase three-site branch of the form (at 
the AC-limit) (0,..., 0, — ~ \[^)^ 0,... 0). Here the stability matrix is the same as in the in-phase case. 

























but with an opposite sign, hence the eigenvalue predictions of the theory of Eq. (USD are the same as in the former case, 
but along the real axis, as opposed to along the imaginary one. As we have seen multiple times with real eigenvalues, 
the predictions are fairly accurate for small e, but for large values of the parameter, higher orders take over and lead 
the pairs to return to the origin. Here, both pairs return to the origin around e = 0.095, the point of the termination 
of the branch. This end point is intriguingly the same as the termination point of both (0,..., 0, W 0,... 0) 


and (the two mirror image installments) (0,..., 0, y ~\J ... 0) and (0,..., u, — y 

effect, we see a rather unusual bifurcation scenario here, which appears as a sort of “double pitchfork”. Namely, 
there are two pairs of eigenvalues involved (hence the “double” designation). For the branch with a single excited 
site, these eigenvalues both come from the imaginary side (bifurcating from the continuous spectrum), while for the 
three-site out-of-phase branch, they both come from the side of the real axis. For the asymmetric branches, the two 
pairs are split with one on the real and one on the imaginary axis. Hence, the bifurcation effectively involves a highly 
symmetric pair of subcritical pitchforks, ultimately leading to the termination of all 4 associated branches. 


, 0 ,- 
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FIG. 5: 


( 0 ,..., 0 ,- 

e = 0.095. 


Same 


9(-l)’ 


for all 



the 

.. 0 ) 


previous branches, but now for the configuration of the form 
at the AC limit, which bears two real eigenvalue pairs and terminates at 


We now briefly explore the “extended” branch in which at e = 0, all t he sites are excited. Given the algebraically 
growing structure of the nonlinear prefactor, the form Vn = \/TJgJn) provides a decaying wave profile. Firstly, 
it is interesting to note here that for all the values of e considered in our computation, this profile was found to 
persist, suggesting, similarly to [s^, that this solution may persist all the way to the continuum limit. It is generally 
worthwhile to iterate here that we found that configurations with progressively larger support were found to persist 
for larger intervals of e values. This is entirely contrary to what is known e.g. for the standard homogeneous focusing 
case (see [4l| for a relevant discussion), where the more localized configurations are the ones eventually persisting 
all the way to the limit, while all others disappear through suitable bifurcations. For the relevant extended solution 
presented in Fig. [6l it is worthwhile to also touch upon its spectrum. Interestingly, since all the sites are excited at 
e = 0, each of them is also associated with a zero pair. Hence, all eigenvalues are initially at the origin and bifurcate 
from there. The result is the apparent discrete spectrum in the right panel of Fig. [6l whereby the eigenvalue pairs 
parabolically grow as e increases. The detailed stability properties of such a configuration merit separate investigation, 
but suffice it to mention for present purposes that the configuration was found to be stable for all the considered 
values of the coupling strength. 

A summary of the different types of states that are examined above is provided in Table IH together with the 
bifurcations leading to their termination and the associated (approximate) critical points. Additionally, in Fig. [71 
we offer an alternative diagnostic that can also be meaningfully used to detect the relevant bifurcations and branch 
collisions. In particular, we show J\f = IwnP as a function of e, which allows to monitor the continuation of 
the solutions for different values of coupling strength parameter. The left panel shows the branches A, C, and E of 
the table, namely the single-site, two out-of-phase and three adjacent out-of-phase sites which collide in the double 
pitchfork bifurcation around e = 0.095, while the right panel illustrates branches B and E, namely the in-phase two- 
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FIG. 6: The profile (left panel) and the stability (right panel) again for e = 0.08 but now for the “extended” configuration in 
which all sites are excited at e = 0 according to: The configuration is generically stable as is also illustrated by 

the spectrum and its zoom-in inset in the right panel. 


site branch and the one with the form (0 ,..., 0, near the e = 0 limit, which, in turn, 

collide and disappear in a saddle-center bifurcation around e = 0.091. 




FIG. 7: The left panel shows the collision of the single-site branch (0,..., 0, ^sM’ (blue solid line) with the 

two out-of-phase site branch (0,..., 0, 0,... 0) (red dashed line) and the three out-of-phase site branch 

(0,. .., 0, gC^TJ’ \l~\llU)' 0,... 0) (green dash-dotted line), through monitoring the dependence of their respective 
powers M = as a function of e. The right panel is similar but now for the collision of (0,..., 0, \llU)' 0,... 0) 

(blue solid line) with (0,... , 0, 0,... 0) (red dashed line). 


Finally, we now turn to the dynamical exploration of the evolution of the unstable configurations in the space-time 
numerical experiments of Figs. 1^1^ and ITOl These are all performed for the case of e = 0.08 used previously to showcase 
the solution profiles. Given the similarity of the profiles of the different branches (and the bifurcations elucidated 
above), we only show three out of the hve few-site excited branches (recall that the single-site excited branch, as well 
as the extended profile branch are stable throughout their respective regimes of existence). Fig. [8] illustrates the case 
of the two-site in-phase excitation branch, Fig. [9] corresponds to the out-of-phase two-site excitation, while Fig. [Till is 
associated with the three-site in-phase excitation. Recall that the mixed phase three-site excitation is rather similar 
in profile to the two-site in-phase, as is the three-site out-of-phase to the two-site out-of-phase for this value of e. 

In all three cases, the two panels, respectively, demonstrate the space-time evolution of the contour of the solution 
magnitude and its difference (again in magnitude) from its initial spatial profile. The former provides a sense of the 
dynamics, while the latter also gives a glimpse of the type of instability that results in it. Interestingly, in all the 
cases we see a rather similar evolution, i.e., over time while the dynamics does not appear to definitively settle to an 
asymptotic state, it does seem to expand its spatial extent, lending further support to the idea that configurations 
with more excited sites are favored in the present setting. On the other hand, we do also detect some differences 
between the different cases. In particular, the oscillatory instabilities of Fig.|5]and Fig. [TUI bear a much weaker growth 
rate (as is typically the case for oscillatory instabilities in comparison to exponential ones), and thus require a far 
longer (by an order of magnitude, which roughly mirrors the corresponding difference in growth rates) time interval 
to manifest themselves in comparison to the rapidly developing exponential growth of Fig. jU) In the former cases, 
the right panel appears to also mirror the oscillatory nature of the instability at its dynamical onset i.e., there is an 
interval of oscillatory growth as is expected by the complex nature of the unstable eigenvalues associated with these 
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cases. It is also relevant to point out that in these cases, to seed the instability a random (uniformly distributed) noise 
has been added to the initial condition, while in the case of Fig. [^this was not necessary (i.e., numerical round-off 
error was rapidly -exponentially- amplified in the latter setting). 



FIG. 8: The left panel shows the space (n)- time (t) evolution of the modulus of the two-site, in-phase solution (of the form 
(0,. .., 0, 0,... 0) at the AC limit). The right panel shows the difference of the magnitude of the solution from 

the magnitude of its corresponding initialization. The coupling strength here is chosen as e = 0.08. 



AC limit). 



IV. CONCLUSIONS & FUTURE CHALLENGES 

In the present work, we have explored a setup of increasing interest in the theory of nonlinear waves in lattices and 
continua, namely the emergence of bright solitary waves in defocusing nonlinear media, in the presence of a spatially 
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inhomogeneous nonlinearity profile. Our specific interest here was to explore the lattice setting, of particular potential 
relevance to waveguide applications. The perspective utilized was that of the so-called anti-continuous limit which 
enabled a systematic theoretical analysis, perturbatively close to that limit. This provided not only a roadmap on the 
available coherent structures, but more importantly a handle on their expected stability properties. The analytical 
results obtained by means of this approach were fully corroborated by detailed numerical existence and spectral 
computations. The latter additionally revealed the bifurcation type scenarios that emerge, as well as provided an 
understanding on which states may be favored in such a setting. We also used a number of proof-of-principle numerical 
simulations in order to explore the dynamical evolution of potentially unstable states. 

We believe that these efforts will provide further insight on the relevant phenomenology and will also give a significant 
amount of motivation for their exploration in experimental setups in nonlinear optics that presently appear to be well 
within reach. Further theoretical efforts could focus on a variety of settings. It would be interesting for example 
to provide an analytical characterization of the spectral operators and the stability of the extended (stable) state 
that we discussed herein, as well as to explore the similarities and differences (existence, stability and dynamics-wise) 
of different “profiles” of the inhomogeneous nonlinearity, such as the exponential one previously studied in vs. 
power-law (as e.g. in the case example considered herein). Our analysis, to the extent possible herein, was kept very 
general, and clearly some features (like the decay of the extended state) will accordingly differ, but if some broad 
qualitative statements could be made along these lines, it would be especially useful (including in designing relevant 
experiments). Lastly, and perhaps most importantly exploring such systems in higher dimensions and identifying the 
impact on such inhomogeneous nonlinearities on different kinds of structures, including vortical ones would be an 
especially relevant theme for future investigations. Efforts in this direction have been recently initiated in 2d settings 
(see e.g. 0) and may well be relevant to extend also to 3d case examples. 
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